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Self-organized criticality [T| is one of the key concepts to describe the emergence of com- 
plexity in natural systems. The concept asserts that a system self-organizes into a critical 
state where system observables are distributed according to a power-law. Prominent ex- 
amples of self-organized critical dynamics include, e.g., piling of granular media [2], plate 
tectonics |3j and stick— slip motion [4j. Critical behavior has been shown to bring about op- 
timal computational capabilities [5], optimal transmission |6| and storage of information ff). 
and sensitivity to sensory stimuli (Sl [Ol llOj . In neuronal systems the existence of critical 
avalanches was predicted in a paper of one of the present authors and observed exper- 
imentally by Beggs and Plenz [6], I12L I13| . Nevertheless, while in the experiments critical 
avalanches were found generically in the sense of genuine self-organized criticality, in the 
model of Ref. they only show up, if the set of parameters is fine-tuned externally to a 
critical transition state. In the present paper we demonstrate analytically and numerically 
that by assuming (biologically more realistic) dynamical synapses [14j in a spiking neu- 
ral network, the neuronal avalanches turn from an exceptional phenomenon into a typical 
and robust self-organized critical behavior if the total resources of neurotransmitter are 
sufficiently large. 

In multi-electrode recordings on slices of rat cortex and neuronal cultures [6l [12] neuronal avalanche 
were observed whose sizes were distributed according to a power-law with an exponent of -3/2. The 
distribution was stable over a long period of time. Variations of the dynamical behavior are induced by 
application or wash-out of neuromodulators. Qualitatively identical behavior can be reached in models 
like [m [15] by variations of a global connectivity parameter. In these models, criticality only shows up, 
if the interactions are fixed precisely at a specific value or connectivity structure. 

Here we study a model with activity-dependent depressive synapses and show that existence of several 
dynamical regimes can be reconciled with parameter-independent criticality. We find that synaptic de- 
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Figure 1: Distribution of avalanche sizes for different coupling strengths a. At a < 1.3 small avalanches 
are preferred yielding a subcritical distribution. The range of connectivity parameters near a = 1.4 
appears critical. For a > 1.6 the distribution is supercritical, i.e. a substantial fraction of firing events 
spreads through the whole system. These results are shown for TV = 300, v = 10, u ~ 0.2, P^'^ = 0.025, 
other parameter values are discussed below. 

pression causes the mean synaptic strengths to approach a critical value for a certain range of interaction 
parameters, while outside this range other dynamical behaviors are prevalent, cf. Fig.[TJ We analytically 
derive an expression for the average coupling strengths among neurons and the average inter-spike inter- 
vals in a mean-field approach. The mean field approximation is applicable here even in the critical state, 
because the quantities that are averaged do not exhibit power-laws, but unimodal distributions. These 
mean values obey a self-consistency equation which allows us to identify the mechanism that drives the 
dynamics of the system towards the critical regime. Moreover, the critical regime induced by the synaptic 
dynamics is robust to parameter changes. 

Consider a network of N integrate-and-fire neurons. Each neuron is characterized by a membrane po- 
tential Q < hi{t) < 9. The neurons receive external inputs by a random process S {1, . . . , N} which 
selects a neuron {t) = i at a rate t and advances the membrane potential hi by an amount P"^^ . Each 
neuron integrates inputs until it reaches a threshold 9. As soon as hi[t) > 9 the neuron emits a spike 
which is delivered to all postsynaptic neurons at a fixed delay r^. The membrane potential is reset by 
^ii^tp) ~ ^i(^sp) ~ ^- For simplicity we will assume in the following that 9 = 1. Super-threshold activity 
is communicated along neural connections of a strength proportional to Jij to other neurons and may 
cause them to fire. In this way an avalanche of neural activity of size L > 1 is triggered. More precisely, 
an avalanche is a period of activity that is initiated by an external input and that terminates when no 
further neuron becomes activated. We define the size of the avalanche as the number of participating 
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neurons. The dynamics of the membrane potential is described by the following equation 



(1) 



In studies of self-organized criticality typically a separation of time scales is assumed which enters in 
Eq-fflby the condition ^ r. It allows us to assume that external input is absent during avalanches. 
Later, in the discrete version of the model, r will play the role of the temporal step size. The variables 
Jij are subject to the dynamics 



which describes the amount of available neurotransmitter in the corresponding synapse [T^ . Namely, if 
a spike arrives at the synapse, the available transmitter is diminished by a fraction u, i.e. the synaptic 
strength decreases due to the usage of transmitter resources. If the presynaptic neuron is silent then the 
synapse recovers and and its strength Jij approaches the value a/u at a slow time scale tj — tvN with 
1 < V <^ N . Thus, the maximal strength of a connection is determined by the parameter a and can be 
observed only when the synapse is fully recovered. The behavior of the network is determined, however, 
by the averaged synaptic strength which will be denoted by (Jy) with the average taken with respect to 
the distribution of inter-spike intervals. In order to obtain our main result we will calculate this effective 
value and use it in a static network. The uniform strengths of the static network are denoted by ag. 

If the external drive and the synaptic weights are small, the activity of the network consists of short burst- 
like events which are initiated by a particular external input. The firing events are separated by relatively 
long relaxation intervals when external inputs are integrated. We may thus be tempted to assume J ~ ^ 
before any spiking event. In general, however, we must take into account that the efficacy of a synapse 
varies in a usage-dependent way which compensates large levels of network activity. Depending on the 
synaptic strength the network can produce a rich repertoire of behaviors. In Fig. [H we show examples 
of avalanche size distributions for various values of a. For small values of a, subcritical avalanche- 
size distributions are observed. This regime is characterized by a negligible number of avalanches that 
extend to the system size. Near a„ the system has an approximate power-law avalanche distribution for 
avalanche sizes almost up to the system size where an exponential cut-off is observed. The mean-squared 
deviation from an exact power law is shown in Fig. [2l Finally, the distribution of avalanche sizes becomes 
non- monotonous when a is well above the critical value a„- 

In preHminary numerical studies we had assumed a model with facilitating and depressing synapses [17] . 
Here we conclude that facilitating synapses are not necessary to evoke self-organized critical avalanches 
in spiking neural networks, depressing synapses are sufficient. This is in Hne with the observation [18] 
that synapses that connect excitatory neurons are largely depressive. To identify the parameters of 
the avalanche size distribution it is sufficient to determine the average synaptic strength: As seen in 
Fig. [3] both the power-law exponent and the mean-squared deviation from the power-law are the same 
for networks with dynamical synapses and networks with static synapses if the strength of the static 
synapses is chosen as ao = u {Jij)- In order to calculate the average synaptic strength analytically, we 
consider in addition the neural inter-spike intervals A'**'. On the one hand, if the inter-spike intervals are 
short then the synapses have a short time to recover and the average synaptic strength resides at a low 
level. On the other hand, large synaptic strengths lead to long avalanches and to large input to neurons 
during the avalanches, which tends to shorten the inter-spike intervals. 




(2) 
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Figure 2: The range of connectivity parameters which cause a critical dynamics extends with system 
size. The mean-squared deviation from the best-matching power-law is plotted in dependence of the 
connection strengths ag and a for static synapses and depressive synapses, respectively. Blue circles, 
squares and triangles stand for networks with dynamical synapses and system sizes = 1000, 700, and 
500, respectively. Red symbols represent the static model. Note that the minimum of all curves depends 
of the network size [TT]. The inset (same symbols) shows the lengths of the parameter intervals where 
the deviation from the best-matching power-law is smaller than an ad-hoc threshold (A7 = 0.005). 
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Figure 3: Rescaling of depressive synapses. A network with static synapses of uniform strength ao has 
the same statistical properties as a network with dynamic synapses if ao is fixed at the average synaptic 
strength of the dynamical case, i.e. if ao := {J)u. The mean-squared deviation A7 from the best-matching 
power-law is shown as a function of the synaptic strength ao for static synapses (red symbols) and the 
mean synaptic strength for dynamical synapses (blue symbols), respectively. The inset (same symbols) 
shows the exponent 7 of the best-matching power-law in the two cases. Parameters are N ~ 100, ly = 10, 
u = 0.2. 
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Figure 4: The graphical solution of Eqs.[3]estabHshes a functional relation between the average synaptic 
strength and inter-spike interval. It is obtained by the intersections of the solutions of Eq. [12] (dashed 
Hue) and Eq. [7] (solid Hues) for a = 1.3, ... , 2.0 in steps of 0.1 (from right to left). This solution agrees 
well with the results of simulations (circles) of a network with the same values of a. Parameters are 
TV = 500, v = lQ,u = 0.2. 

This trade-off determines the expected values of the synaptic strengths and the inter-spike intervals which 
are reaHzed by the dynamics of the network. In order to express this reasoning more formally, we solve the 
dynamical equations ^ and ([2]) based on a stationarity assumption for both the synaptic strengths and 
the inter-spike interval. Neither of these quantities has a power-law distribution and their first moments 
exist. In Methods we derive expressions of the mean synaptic strength {Jij) and the mean value of the 
inter-spike intervals distribution (A'^'). The stochastic dependency of the two quantities is reflected in 
a mutual dependence of their averages. Each of the dependencies is derived analytically which allows 
us to formulate the self-consistency of the stationarity requirement as the simultaneous solution of the 
mean-field equations 

{J,,) = {J,,) ( (A'^')) and (A-) = (A^') ( (J,,-)) (3) 

which can be determined graphically from the intersections of the solutions of Eqs. [7]and[l2l cf. Fig. IH 
The mean-field solution is confirmed by direct network simulations that are are represented by the circles 
in Fig. m The solution is unique for any a. The stationary distribution is less sensitive to changes of 
the parameter a near the critical value of the synaptic strength than further away from it. This brings 
about the large critical region for the model with depressive synapses, cf. Fig. [2j 

Furthermore we want to discuss the stability of the solution of the self-consistency equation ([3]) . If we 
apply a perturbation A J to all synapses at time tp such that for each i,j Jij = Jij + A J, we can show 
with srane simple computations, that before the next spike the synaptic strengths are on average smaller 
than Jij. In the simulated system the average synaptic strength is driven back to the fixed point by a 
few spikes, such that the solution of ([3]) is indeed stable for the critical state. 

Both the numerical study in Ref. [17] and the analysis presented so far refer to finite systems. In order 
to check whether the trend that is visible in Fig. [2] continues for larger network sizes, we consider the 
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behavior of the mean synaptic strength in the thermodynamic hmit N ^ oo we compute the expectation 
value of the avalanche-size distribution ifTTj) , (L) = jY„(-j^]^-)^^ , and insert it into the self-consistency 
equation ^ 



In the limit iV — > oo we should scale the external input /™* ~ A^~"' and w > 0. We now distinguish 
the following cases, a) If (A"^') ^ iV^+^ and e > w then the right hand side (r.h.s.) of ^ tends to 
— oo, while the l.h.s. is always larger than 0. b) If (A'*'') ~ iV^+'^ and < e < u; then the r.h.s. of ^ 
tends to 1 (or if e = if) while the l.h.s. a, hence a solution is only possible if a = 1 and in this case 
u{Jij) -> 1. c) If (A™) ^ iV^~^ and e < then the r.h.s. of ^ tends to 1, while the l.h.s. approaches 
0. d) If (A'^') ^ TV, we can assume that (A"'') = cN + o{N). From ([4]) for a > 1 we can find the unique 
solution c = —ly (ln(a — 1) — ln(a — 1 + u)). In all cases when the solution exists, u{Jij) 1, which we 
know to be the critical connectivity for the network with statical synapses in the Hmit iV ^ oo. Hence 
in the thermodynamical limit the network with dynamical synapses becomes critical for any a> 1. 

In this paper we have focused on fully connected networks and neurons without leak currents for reasons of 
analytical tractability. We now discuss the results of various generalizations which we have investigated 
numerically. If the network has only partial connectivity, the results stay the same, if the synaptic 
strengths are properly rescaled. In a random network of size N with connectivity probability c, the 
critical parameter a is approximately equal to a'j^/c, where a'§ is obtained from the critical parameter 
region of the fully connected network of size c x N. If the connections in a partially connected random 
network are not chosen independently (e.g. "small-world" connectivity [19]) one finds even more accurate 
power-laws than for the independent case with the same average connectivity. A similar phenomenon 
occurs in the grid network [l6j which has been used to model criticality in EEG recordings. 



If in Eq. [T]we add a leak term, which is present in biologically reaHstic situations 

TV 

N 



1 ^ 

h{t) = ~T-'Ht) + C + <5.,«,(t)/''^' + ""J^o^ - ^sp - ^d) , (5) 

we find numerically that the distribution of the avalanche sizes remains a power law for leak time-constants 
up to Ti « 40ms. In ([5]) we included a constant compensatory synaptic current C which depends on r/ 
and summarizes neuronal self-regulatory mechanisms. In this way the probability of the neuron to stay 
near threshold is conserved and avalanches are triggered in a similar way as in the non-leaky case. The 
resulting power-law exponent is sHghtly smaller than —1.5 and reaches values close to —2 for strong 
leakage in simulations of a network of iV = 500 neurons. 

In summary, we have presented an analytical and numerical study of spiking neural networks with dy- 
namical synapses. Activity-dependent synaptic regulation leads to the self-organization of the network 
towards a critical state. The analysis demonstrates that mean synaptic efficacy hereby plays a crucial 
role. We explained how the critical state depends on the maximally available resources and have shown 
that in the thermodynamic limit the network becomes critical for any a > 1, i-e. that criticality is an 
intrinsic phenomenon produced by the synaptic dynamics. The mean field quantities are in very good 
agreement with simulations and were shown to be robust with respect to perturbations of the model 
parameters. 
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Methods 

In this section we give an explicite derivation of the self-consistency relation ([3]) . Solving Eq. [2] between 
two spikes of neuron j we find 

J,,iq) = ^ (l - (l - ^ J,,(t+)) e-(*-*^)/-) , (6) 

where the synaptic strengths immediately before and after a spike of neuron j at time ts are denoted by 
Jij{tj) and Jij{tf), respectively. Within a short interval containing the spike, decreases by a fraction 
u such that Jij{ti) = (1 — u) Jij{t^). 

The average synaptic strength (Jy) is the expectation value of Jij{tj). Analogously, (A"'') refers to 
the inter-spike interval A"^'. The random variables Jij{t~) and A"^' both depend on the distribution of 
external inputs and are thus related. The self-consistency conditions ^ express this relation for the 
respective averages. Assuming that (Jy ) depends essentially only on the mean inter-spike interval, we 
obtain from 

u 1 — (1 — u)e \^ / 
where the average is performed in discrete time with step width r, i.e. tj = vN . 

We are now going to establish a relation between P (A"") and the inter-avalanche interval distribution 
Q (A'^') . It can be shown that the neuronal membrane potentials before an avalanche are uniformly 
distributed on the interval where is a lower bound of hi{tsp) — with eat for ~^ cxo. 

Under these conditions, (3(A'^') has a geometric distribution 



Q (A'^M = 1 . (8) 

Let kj be the number of avalanches between two spikes of the neuron j. A mean-field approximation 
relates the averages of the distributions of inter-spike and inter-avalanche intervals 

(A'^') = (fc)(A'"'). (9) 

The average inter-avalanche interval is easily computed from ([8]) 

i^'n = (10) 

In order to calculate the average number of avalanches between two spikes of a neuron, we compute 
the time to reach the threshold by accumulating external inputs and spikes from other neurons during 
avalanches, i.e. 
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where (L) is the mean avalanche size and is the probabiUty that neuron j is receiving an input. The 
distribution of avalanche sizes can be computed analytically for a network with static synapses of strength 



In the case of dynamical synapses we apply a mean field approximation and set ao ^ u (Jij) in ifTTj) . This 
allows us to compute (L) as a function of (u {Jij)). 

Combining ([9]), (fTOj) . and ifTTj) we obtain a relation between the inter-spike interval and the average 
synaptic strength. 

/A^s,^0_^ ^ (12) 

The self-consistency equations ^ arise from ([7]) and l|12p . Its solution is obtained by numerical analysis 
of the two independent relations ^ and (fT2l) . cf. Fig.Sl 
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